% created by tyy on 2023/11/07
function [x_out,f_out,cnt,time,code] = f_steffensen_method(equation,range,tolerance,maxIteration)
% Steffensen 法
syms x;
% iteration parameters
a = sum(range)/2; % initial point

% Newton Iterative Method
k_x = zeros(maxIteration,1);
f_x = zeros(maxIteration,1);
tic;
x_0 = a;
for i = 1:maxIteration
    f_x(i,1) = subs(equation,x,x_0);
    k_x(i,1) = (f_x(i,1) - subs(equation,x,x_0-f_x(i,1))) / f_x(i,1);
    delta = - f_x(i,1) / k_x(i,1);
    % delta  = - f_x(i,1)^2 / ( subs(equation,x,x_0-f_x(i,1)) - f_x(i,1) );
    x_0 = x_0 + delta;
    % x_new = x_value - subs(equation,x,x_value)/k_x(i,1);
    if abs(double(subs(equation, x, x_0))) < tolerance
        code = 0;
        break;
    end
end
time = toc;

if (i>=maxIteration) 
    code = -1; 
end
x_out = x_0;
f_out = subs(equation,x,x_out);
cnt = i;